% Function fzeroall finds all zeros of function f
% in range [span(1), span(2)]
% No more than 1 zero per n-th part of span

function z = fzeroall(f, span, n, varargin)
    
    if (nargin == 4) % optional vector_acceptable is present
        vector_acceptable = varargin{1, 1};
    else
        vector_acceptable = true;
    end

    a = span(1);
    b = span(2);
    step = (b - a) / n;
    
    x = linspace(a, b, n+1)';
    
    if (vector_acceptable)
        % f accepts matrix argument
        y = f(x);
    else
        % f does not allow matrix argument
        y = zeros(size(x));
        for i = 1:size(x, 1)
            y(i) = f(x(i));
        end
    end
    
    s = (y(1:n) .* y(2:n+1)) <= 0;
    
    indexes = find(s > 0.5);
    z = zeros(size(indexes));
    
    for i = 1:size(indexes, 1)
        idx = indexes(i);
        
        m = fzero(f, a+step*[idx-1, idx]);
        
        z(i) = m;
    end

end